scRNA-seq

## load required packages
library(Seurat)
library(cowplot)
library(dplyr)
library(ggplot2)
library(DT)
library(paletteer)
library(forcats)


orig.ident_fils = paletteer::scale_fill_paletteer_d("ggsci::nrc_npg") 

DEG Analysis overview

Differentially Expressed genes between TRT and CTL

Number of Cells

res = "sample"
obj.tmp@meta.data %>% ggplot(aes(!!sym(res), fill=!!sym(res))) + 
  geom_bar(alpha=0.7, color="grey5", size=0.1) +
  geom_text(stat="count", aes(label= ..count..), vjust=-0.5, size=3) +
  orig.ident_fils + 
  xlab("") + 
  theme_classic() +
  theme(legend.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust=1)) 

Analysis

Volcano plot

# Define funtion 
#id1 = "treatment"
#id2 = "control"
#logfc = 0
deg.two.groups = function(obj.srt = obj.srt, treatment, control,logfc=0){
  Idents(obj.srt) = 'sample'
  markers <-FindMarkers(
    obj.srt,
    logfc.threshold = logfc,
    ident.1 = treatment,
    ident.2 = control, slot= "data")
}
deg = deg.two.groups(obj.srt = obj.tmp,
                    treatment = "TRT", 
                    control = "CTL")
t= paste0(paste0("TRT"), " / ",paste0("CTL"))
deg %>% ggplot(aes(avg_log2FC, -log10(p_val))) + 
  geom_point(size=0.1) +
  geom_vline(xintercept = 0, size=0.1) +
  theme_bw() + ggtitle(t)

Volcano plot with DEG information

deg = deg %>% mutate(DE= ifelse(p_val <= 0.05 & avg_log2FC >= log2(1.2), "UP",
                                ifelse(p_val <= 0.05 & avg_log2FC <= -log2(1.2), "DN",
                                "no significant")))
deg$DE = factor(deg$DE, levels = c("UP","DN","no significant"))
up = nrow(deg[deg$DE == "UP", ])
dn = nrow(deg[deg$DE == "DN", ])
deg %>% ggplot(aes(avg_log2FC, -log10(p_val), color=DE)) + 
  geom_point(size=0.5, shape=19, alpha=0.7) +
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), size=0.1, color="grey88") +
  geom_hline(yintercept = -log10(0.05), size=0.1) +
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_bw() +
  annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up), 
           hjust = 1.1, vjust = 2, size = 5, color = "red") +
  annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn), 
           hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
  theme_bw() + ggtitle(t)

significance : p value <= 0.05
log2FC >= 1.2

DEG table

deg %>% DT::datatable()
# GSEA related functions 
library(clusterProfiler)

perform_GSEA <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    # Check the name of log2fc related 
    if ("avg_log2FC" %in% names(res)) {
      df <- res$avg_log2FC
    } else if ("log2FoldChange" %in% names(res)) {
      df <- res$log2FoldChange
    } else {
      stop("Neither avg_log2FC nor log2FoldChange columns found in the data frame.")
    }
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  result <- x@result %>% arrange(desc(NES))
  result <- result[, c('NES', 'pvalue', 'p.adjust', 'core_enrichment', 'ID')]
  return(result)
}

# GSEA Plot 
gsea_nes_plot <- function(gsea.res, title, color="pvalue") {
  gsea.res = gsea.res %>% mutate(sig=ifelse(pvalue <= 0.05, "p value <= 0.05", "p value > 0.05"))
  # basic plot
  p <- gsea.res %>%
    ggplot(aes(reorder(ID, NES), NES)) +
    geom_col(aes(fill=!!sym(color)), color="grey5", size=0.15, alpha=0.8) +
    coord_flip() +
    labs(x="Pathway", y="Normalized Enrichment Score", title="GSEA") +
    theme_classic() +
    theme(axis.text.x = element_text(size=5, face = 'bold'),
          axis.text.y = element_text(size=6, face = 'bold'),
          axis.title = element_text(size=10)) +
    ggtitle(title)
  
  # color by color input type
  if (color == "pvalue") {
    p <- p + scale_fill_gradient(low = 'orangered', high = '#E5E7E9')
  } else if (color == "sig") {
    p <- p + scale_fill_manual(values = c("orangered", "#E5E7E9"))
  }
  return(p)
}
# Pathway selection : HALLMARK 
hallmark <- msigdbr::msigdbr(species = "Homo sapiens", category = "H") %>% 
  dplyr::select(gs_name, gene_symbol) 
gsea.res = perform_GSEA(res = deg, ref = hallmark, pvalueCutoff = 1)

GSEA plot colored by p-value

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "pvalue") 

GSEA plot colored by significance

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "sig")

GSEA information table

gsea.res %>% DT::datatable()